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ABSTRACT 

This  report  describes  an  implementation  of  a  Godunov-type  solver  for  gas  dy¬ 
namics  equations  in  MATLAB®.  The  main  attention  is  paid  to  providing  a 
generic  code  that  can  be  easily  adapted  to  particular  problems  in  one,  two  or 
three  dimensions.  This  is  achieved  by  employing  a  cell  connectivity  matrix 
thus  allowing  one  to  use  various  structured  and  unstructured  meshes  with¬ 
out  modification  of  the  core  solver.  The  code  has  been  thoroughly  tested  for 
MATLAB  Version  7.6  (Release  2008a). 
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Solving  multi-dimensional  problems  of  gas  dynamics  using 

MATLAB® 

Executive  Summary 

In  many  circumstances  it  is  required  to  simulate  blast  propagation  in  complex  three- 
dimensional  domains  in  order  to  estimate  pressure  and  temperature  fields  at  some  distance 
from  the  source  of  explosion.  Though  this  is  the  classical  problem  of  gas  dynamics  whose 
solution  is  implemented  in  many  commercial  software  packages,  the  use  of  such  pack¬ 
ages  is  sometimes  limited  to  a  particular  platform,  specific  mesh  type  and  proprietary 
input /output  file  formats,  and  requires  a  long  learning  curve  associated  with  input  data 
preparation  (pre-processing)  and  visualisation  of  the  obtained  results  (post-processing). 
At  the  same  time  commercial  software  may  not  be  open  to  model  extension  and  may  not 
be  easily  embedded  into  other  systems. 

This  report  describes  an  implementation  of  a  Godunov-type  solver  for  gas  dynamics 
equations  in  MATLAB®.  The  main  attention  is  paid  to  providing  a  generic  code  that 
can  be  easily  adapted  to  particular  problems  in  one,  two  or  three  dimensions.  This  is 
achieved  by  employing  a  cell  connectivity  matrix  thus  allowing  one  to  use  various  struc¬ 
tured  and  unstructured  meshes  without  modification  of  the  core  solver.  The  code  has  been 
thoroughly  tested  for  MATLAB  Version  7.6  (Release  2008a). 

The  work  has  been  performed  within  the  “Blast  Modelling  for  Cordon  Assessment  and 
Bomb  Scene  Examination”  research  project  NS  07/002. 
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1  Introduction 


In  many  circumstances  it  is  required  to  simulate  blast  propagation  in  complex  three- 
dimensional  domains  in  order  to  estimate  pressure  and  temperature  fields  at  some  distance 
from  the  source  of  explosion.  Though  this  is  the  classical  problem  of  gas  dynamics  whose 
solution  is  implemented  in  many  commercial  software  packages,  the  use  of  such  pack¬ 
ages  is  sometimes  limited  to  a  particular  platform,  specific  mesh  type  and  proprietary 
input /output  file  formats,  and  requires  a  long  learning  curve  associated  with  input  data 
preparation  (pre-processing)  and  visualisation  of  the  obtained  results  (post-processing). 
At  the  same  time  commercial  software  may  not  be  open  to  model  extension  and  may  not 
be  easily  embedded  into  other  systems. 

MATLAB®  provides  a  nice  environment  for  numerical  simulation  of  real-world  prob¬ 
lems  with  integrated  visualisation,  powerful  scripting  framework,  and  fast  algorithms  im¬ 
plemented.  Many  books  are  devoted  to  the  application  of  MATLAB  to  various  problems 
of  science  and  engineering  (see  e.g.  [Danaila  et  al.  2007,  Klee  2007,  Gilat  2008]). 

In  this  report,  two  compact  MATLAB  files  (M-files)  constituting  the  development  of  a 
Gas  Dynamics  Toolbox  (MATLAB  library)  are  presented  which  are  sufficient  to  run  from 
a  user-defined  script  and  simulate  gas  dynamics  problems.  A  robust  Godunov- type  solver 
[Godunov  1959,  Richtmyer  &  Morton  1967,  Holt  1977,  Toro  1999]  based  on  the  exact 
Riemann  solver  for  flux  calculation  between  control  cells  is  employed.  Though  there  exist 
many  advanced  schemes  for  simulation  of  gas  dynamics  equations  with  higher  accuracy  (see 
e.g.  [van  Leer  1979]),  the  first-order  accurate  Godunov  solver  provides  a  highly  desirable 
monotone  numerical  scheme  and  is  easily  extendable  to  multi-dimensional  problems.  The 
history  of  discovery  of  this  algorithm  based  on  deep  insight  into  the  physics  of  shock  and 
rarefaction  waves  is  given  in  [Godunov  1999]. 


2  Elements  of  gas  dynamics 


Let  w  be  a  control  volume  with  piecewise  smooth  boundary  du)  oriented  by  inward 
normal  unit  vector  n.  The  classical  conservation  laws  of  mass,  momentum  and  energy  for 
gas  flow,  written  for  arbitrary  w,  are  as  follows 

A  J  p^v  =  J  •  ndA,  (1) 

dm 


dt 


I 


LJ 


pvdK 


/  (yov  V  •  n -|-pn)  dA  , 
duo 


dt 


LJ 


+  P 


V  •  n  dA  . 


(2) 

(3) 


Here  t  is  time,  p  density,  v  velocity,  p  pressure,  and  e  specific  internal  energy.  Symbols  dV 
and  dA  denote  volume  and  area  elements,  respectively.  The  effect  of  gravity  is  ignored. 
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Since  the  control  volume  u)  is  arbitrary,  the  integral  conservation  laws  (l)-(3)  imply 
the  following  differential  conservation  equations 


|  +  V.(pv)  =  0, 


(4) 


*9(/9v) 

dt 


+  V  •  (pv  (g)  V  +pG)  =  0  , 


Wt 


p 


+  V- 


p 


+  p 


=  0. 


(5) 

(6) 


where  V  denotes  the  gradient  operator,  (g)  the  tensor  product,  and  G  the  metric  tensor. 
The  conservative  form  of  gas  dynamics  equations  admits  discontinuous  solutions  such 
as  shock  waves  and  contact  discontinuities.  In  this  case  temporal  and  spatial  derivatives 
have  to  be  considered  in  the  generalised  sense.  For  continuous  solutions  such  as  rarefaction 
waves  the  conservative  Euler  equations  (4)-(6)  can  be  further  simplified  to 


dp 


+  V  •  Vp  +  p  V  •  V  =  0  , 


(7) 


d'v 

p(— +  v-Vv)+Vp  =  0, 


+  v-Ve  +pV-v  =  0. 


de 
dt 

As  a  consequence  of  equations  (7),  (9)  and  the  second  law  of  thermodynamics 


(8) 

(9) 


0ds 


de  +pd 


(10) 


where  6  and  s  are  the  absolute  temperature  and  specific  entropy,  respectively,  the  sub¬ 
stantial  derivative  of  entropy  vanishes: 


ds 

dt 


ds 


+  V  •  Vs  =  0  . 


(11) 


In  other  words,  entropy  at  a  particle  remains  constant  in  a  continuous  gas  flow.  Note  that 
entropy  at  a  particle  undergoing  a  shock  wave  always  increases. 

The  fundamental  thermodynamic  relations  for  ideal  polytropic  gas  result  in  the  fol¬ 
lowing  state  equation  [Courant  &  Friedrichs  1948] 


p  =  (y-  l)pe 


(12) 


where  7  is  the  ratio  of  the  specific  heat  at  constant  pressure  to  the  specific  heat  at  constant 
volume.  Indeed,  by  the  definition,  gas  is  ideal  (thermally  perfect)  if  its  thermodynamic 
state  satisfies  Clapeyron’s  equation 


p  =  RpO 


(13) 
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where  R  is  the  specific  gas  constant  defined  as  the  ratio  of  the  universal  gas  constant 
to  the  molar  mass  of  gas  (or,  equivalently,  the  ratio  of  the  Boltzmann  constant  to  the 
average  mass  of  the  gas  molecule).  From  the  second  law  of  thermodynamics  (10)  it  is 
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straightforward  to  deduce  that  the  specific  internal  energy  of  ideal  gas  is  a  function  of 
the  absolute  temperature  alone.  Ideal  gas  is  polytropic  (calorically  perfect)  if  e  is  a  linear 
function  of  9,  namely 

e  =  Cv9  (14) 

where  c„  is  the  specific  heat  at  constant  volume.  In  this  case  the  specific  enthalpy  h, 
defined  by  the  expression  h  =  e+p/ p,  is  also  a  linear  function  of  the  absolute  temperature 

h  =  CpO  (15) 

where  Cp  =  Cv  +  R  is  the  specific  heat  at  constant  pressure.  These  formulae  result  in 
expression  (12)  with  7  =  Cpjcv  Note  that  the  sound  speed  a  of  a  polytropic  gas  is  given 
by  the  formula 


Let  us  introduce  the  volumetric  density  of  momentum,  u  =  pv,  and  the  volumetric 
density  of  total  energy,  £  =  p  .  In  terms  of  these  conserved  variables  the 

conservation  laws  (l)-(3)  take  the  succinct  form 


r,,e]dV  =  l 

^  du) 


1  e  +  p 

u  n,  -  uu  n+pn,  - u  •  n 

P  P 


dyl 


where 


p  =  (7-  1)  e- 


u 


2p 


(17) 


(18) 


according  to  the  equation  of  state  (12). 

Most  numerical  algorithms,  including  Godunov-type  solvers,  are  based  on  the  above 
integral  equations.  The  main  problem  is  related  to  the  correct  determination  of  the  fluxes 
of  mass,  momentum  and  energy  through  the  interface  separating  two  control  volumes  in 
terms  of  the  cell  values  of  the  conserved  variables. 


3  The  Riemann  problem 

The  one-dimensional  Riemann  problem  is  important  for  understanding  the  physics  of 
compression  and  expansion  waves,  and  provides  an  exact  expression  for  fluxes  between  two 
adjacent  uniform  states  of  gas  separated  by  an  infinite  plane.  The  solution  of  the  Riemann 
problem  is  described  in  many  texts  on  gas  dynamics  (see  e.g.  [Courant  &  Friedrichs  1948]). 
In  this  section  the  solution  is  outlined  without  providing  full  details. 

Consider  the  propagation  of  a  shock  wave  or  rarefaction  fan  into  a  uniform  gas  in 
the  direction  of  Cartesian  coordinate  x  with  front  initially  positioned  at  x  =  0.  It  is 
assumed  that  the  gas  flow  does  not  depend  on  the  remaining  Cartesian  coordinates,  and 
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velocity  field  has  the  only  nonzero  component  v  in  the  x-direction.  Denote  given  pre-state 
of  gas  initially  occupying  the  half-space  {x  >  0}  by  {po-,  Po,  Vq)-  Let  us  determine  all  the 
possible  post-states  (p*,  p*,  u*)  of  gas  initially  occupying  the  complementary  half-space 
{x  <  0}  and  connected  with  the  given  pre-state  by  a  shock  wave  or  simple  rarefaction  fan, 
both  propagating  in  the  positive  direction  of  x  relative  to  Vo-  It  turns  out  that  the  post¬ 
states  cannot  be  arbitrary  but  form  a  one-dimensional  manifold  (curve).  This  curve  can 
be  parametrised  by  the  post-state  pressure  p*,  and  therefore  the  density  p*  and  velocity 
u*  are  known  functions  of  p*  and  the  pre-state  {po,  Po,  Vq)-  These  functions  are  piecewise 
analytic,  described  by  two  different  expressions  for  p*  >  po  (Hugoniot  curve)  and  p*  <  po 
(Poisson  curve).  The  Hugoniot  and  Poisson  curves  join  smoothly  at  p*  =  po  up  to  the 
continuous  second-order  derivative  and  form  the  combined  Hugoniot-Poisson  adiabatic 
curve  given  by  the  formulae 


—  =  R 

Po 


7 


(19) 


U*  -  Vo 


=  K 


(Xn 


Here 


R^{q)  =  I 


a+^ 

y  ^  7+1 


7-1 

r+1 

1/7 


V^{q)  =  { 


q-1 


(<?>!) 

(0  <  g  <  1) 

(<?>!) 


(0  <  (7  <  1) 


k  7  —  1 


(20) 


(21) 


(22) 


and  Oo  is  the  pre-state  sound  speed  of  gas.  The  graphs  of  the  functions  R-y{q)  and  V^{q) 
are  displayed  in  Figure  1.  In  acoustic  approximation  of  small  p*  —  Po  (compared  to  po) 
equations  (19)  and  (20)  reduce  to 


P*  -Po 

P*  Po  —  9  ? 


P*  -Po 
po  ^o 


(23) 


The  product  p  a  is  called  impedance. 

Note  that  the  Hugoniot  curve  {q  >  1)  is  derived  from  the  Rankine-Hugoniot  conditions 
at  the  shock  wave,  whereas  the  Poisson  curve  (0  <  (7  <  1)  is  obtained  from  the  exact 
solution  of  gas  dynamics  equations  written  in  terms  of  Riemann  invariants.  The  gas  state 
is  obtained  as  an  explicit  similarity  solution  (depending  on  =  x/t)  to  gas  dynamics 
equations.  For  p*  >  po  the  compression  Riemann  wave  is  given  by  the  piecewise  constant 
profile  of  density,  pressure,  and  velocity 


{P,  P,  v) 


{Po,  Po,  Vo)  {vs  <  C) 

{p*,p*,vX)  {i<Vs) 


(24) 
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(i)  Plots  of  R-y{q). 


(ii)  Plots  of  V^{q). 


Figure  1:  Dimensionless  Hugoniot-Poisson  eurves  for  various  7. 


where 


Vs 


Po  Vo 
po 


(25) 


is  the  shock  speed.  For  0  <  p*  <  Po  the  expansion  Riemann  wave  is  described  by  the 
following  continuous  functions 


(P,  P,  v) 


{Po,  Po,  Vo) 

'  {p{^),p{^),v{0) 

,  {p*,  P*,  V^,) 


{Vo  +  flo  <  0 

{v^  +  <  Vo  +  Oo) 

if,  <V^  +  a*) 


(26) 


where  a*  is  the  post-state  sound  speed  of  gas,  and 


KO  =  Po 


“(0 


2 

7-1 


(27) 


p(0  =  Po 


m 


Po 

v{0  =  ^  -  «(0  > 

T  —  ^ 

«(0  =  - —  {vo  +  ao-C)  . 

7  +  1 


(28) 

(29) 

(30) 


It  is  straightforward  to  check  that,  since  the  state  (p*,  p*,  u*)  belongs  to  the  Poisson  curve, 
the  gas  properties  in  the  expansion  wave,  given  by  (26),  vary  continuously. 


Note  that  equation  (19)  expressed  in  terms  of  pressure  p  as  a  function  of  specific  volume 
1/p  is  actually  called  the  Hugoniot-Poisson  curve.  Conceptually,  the  latter  variables  are 
just  specialised  coordinates  on  the  Hugoniot-Poisson  manifold  of  all  possible  post-states. 


Now  consider  two  arbitrary  uniform  states  of  gas  initially  separated  by  some  plane 
without  assumption  that  the  velocity  vectors  are  perpendicular  to  the  plane.  To  deter¬ 
mine  the  breakdown  of  this  configuration  when  time  increases  the  above  solutions  can  be 
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used.  Indeed,  consider  either  of  half-planes  and  make  Galilean  transformation  to  a  new 
coordinate  system  in  which  the  uniform  velocity  field  is  perpendicular  to  the  plane.  Then 
apply  the  above  analytic  formulae  to  describe  all  possible  post-states  parametrised  by  the 
post-state  pressure.  The  two  post-states  defining  two  Riemann  waves  have  to  be  connected 
to  provide  global  solution.  Since  the  gas  density  and  tangential  velocity  can  be  discon¬ 
tinuous  (due  to  potentially  different  Galilean  transformations),  the  only  kind  of  interface 
between  the  uniform  post-states  can  be  a  contact  discontinuity.  Gontact  interfaces  are 
defined  by  the  conditions  of  continuity  of  pressure  and  the  normal  component  of  velocity. 
These  conditions  are  readily  achieved  by  the  above  explicit  solutions.  Indeed,  the  only 
condition  to  satisfy  is  the  continuity  of  the  normal  component  of  the  post-state  velocity 
as  a  function  of  common  post-state  pressure  p*.  This  results  in  the  following  algebraic 
equation  to  solve 


f{p*)  =  v+ 


V-  +  aj^V^ 


“t"  Cl—  Gy 


=  0 


(31) 


where  the  subscripts  ±  supply  the  two  sets  of  the  pre-state  variables  of  gas  initially  oc¬ 
cupying  half-spaces  x  >  0  and  x  <  0  respectively  {v±  are  the  normal  components  of  the 
pre-state  velocities).  Since  the  combined  Hugoniot-Poission  adiabatic  curve  is  globally 
smooth,  the  fast  converging  Newton-Raphson  method  can  be  applied  to  solve  the  non¬ 
linear  equation  (31).  The  acoustic  approximation  is  frequently  used  as  an  initial  guess  in 
the  iterative  procedure.  Since  /  (p*)  increases  monotonically  and  /(-|-oo)  =  -|-oo,  a  unique 
solution  p*  exists  if  and  only  if  /(O)  <  0  which  is  equivalent  to  the  condition 

2  (a+  -I-  a_) 

V+-V-  <  - ^ - 

7-1 


Some  care  should  be  taken  in  the  situation  of  strong  rarefaction  waves  (/(O)  >  0)  leading 
to  vacuum  condition  [Toro  1999] . 

An  alternative  Riemann  solver  can  be  obtained  by  inverting  equation  (20)  and  solving 
a  non-linear  equation  for  the  post-state  velocity  u* .  The  comparative  performance  analysis 
of  various  exact  Riemann  solvers  is  given  in  [Gottlieb  &  Groth  1988]. 


4  Numerical  model 


Let  {tOi}  (i  =  1, . . . ,  Nc)  constitute  tessellation  of  the  computational  domain  D  where 
Nc  is  the  total  number  of  cells  (control  volumes).  The  state  of  the  gas  dynamics  flow  is 
approximated  by  the  cell  average  values 

Si  =  J  [p,  u,  e]  dG  (33) 

UJi 


and  the  fluxes  by  the  surface  integrals 


duii 


1  ,  £+P 

u  n,  uu  n-|-pn,  - u  •  n 

P  P 


dA 


(34) 
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where 

=  J  dV  (35) 

CJj 

is  the  cell  volume.  According  to  the  conservation  laws  (17)  and  the  equation  of  state  (18), 
the  following  system  of  ordinary  differential  equations  results 

^  =  (36) 

In  order  to  complete  this  system  of  dynamic  equations  one  needs  to  express  the  rates  Ri 
in  terms  of  the  states  Si.  This  is  done  by  employing  the  Riemann  solver  which  allows 
one  to  explicitly  determine  the  fluxes  between  two  cells  in  contact,  provided  that  the  gas 
states  in  the  cells  are  approximated  by  constant  values  and  the  cells  are  separated  by  a 
planar  patch.  In  this  case  the  cell  average  values  Si  are  approximated  by  the  first  order 
quadrature  formulae 

Si  =  [pi,  Uj,  £i]  .  (37) 

Let  us  assume  that  each  cell  ooi  is  a  polyhedron.  The  collection  {cji}  (i  =  1, . . . ,  Nc) 
can  be  potentially  composed  of  polyhedra  of  different  type,  which  do  not  necessarily  share 
facets.  Let  {i^k}  {k  =  be  the  collection  of  all  distinct  atomic  facets  of  the 

polyhedra  where  Nf  denotes  the  total  number  of  facets.  The  term  ‘atomic  facet’  means 
the  following.  If  two  polyhedra  in  contact  do  not  share  a  facet,  only  the  smaller  (atomic) 
facet  is  added  to  the  collection  of  facets  (pk-  For  example,  when  two  tetrahedra  are  attached 
to  one  face  of  a  cube  (regular  hexahedron) ,  only  the  two  triangular  faces  of  the  tetrahedra 
covering  the  square  face  of  the  cube  are  accounted  for.  In  this  case  the  square  face  is  the 
union  of  the  two  atomic  triangular  faces. 

In  order  to  uniformly  treat  boundary  conditions,  it  is  helpful  to  attach  a  ghost  cell  to 
each  boundary  facet.  The  state  of  the  ghost  cell  must  be  specified  by  the  user  to  model 
appropriate  boundary  conditions  [Oran  &  Boris  1987].  This  approach  slightly  increases 
the  total  number  Nc  of  cells  but,  as  a  compensation,  all  the  facets  become  internal  that 
greatly  simplifies  the  logic  of  computation.  It  is  worthwhile  noting  that  periodic  boundary 
conditions  can  be  handled  by  the  cell  connectivity  matrix  alone  without  a  need  to  introduce 
ghost  cells. 

Denote  the  area  of  facet  pk  by  select  a  unit  normal  vector  n^,  and  introduce  &  Nf- 
by-2  facet-to-cell  connectivity  matrix  a  with  entries  defined  as  follows.  If  facet  pk  belongs 
to  the  intersection  n  doji^  and  points  from  wq  to  wq,  then  set  a{k,l)  =  ii  and 
a{k,  2)  =  12-  Mathematically,  the  cell  connectivity  matrix  a  defines  a  directed  graph  (V,  T) 
where  vertices  V  are  cells  (including  ghost  cells)  and  edges  S  are  facets.  In  particular,  a 
basic  theorem  of  graph  theory  [Diestel  2005]  immediately  provides  the  following  relation 

2  if;]  =  ^  deg  (u)  (38) 

v&V 

where  deg  (u)  is  the  degree  (or  valency)  of  a  vertex  v  defined  as  the  number  of  connected 
neighbours,  and  the  standard  symbol  j5j  denotes  the  cardinality  of  a  set  S  (the  number 
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of  elements).  This  relation  can  be  used  to  relate  the  number  of  facets  with  the  total 
number  Nc  of  cells  (including  ghost  cells)  and  the  number  Ng  of  ghost  cells.  For  example, 
if  each  cell  cu*  is  a  n-hedron,  then 

2Nf  =  n  {Nc  —  Ng)  +  Ng  <  n  Nc  ■  (39) 

Indeed,  the  degree  of  each  ordinary  cell  is  equal  to  n  (the  number  of  faces)  and  the  degree 
of  each  ghost  cell  is  equal  to  1  by  construction. 

The  numerical  algorithm  is  straightforward.  First  update  the  states  of  the  ghost  cells 
according  to  the  imposed  boundary  conditions.  For  example,  for  an  ideally  reflective  wall 
or  symmetry  plane,  set  the  ghost  cell  state  equal  to  that  of  the  neighbouring  boundary 
cell  except  for  the  normal  component  of  momentum  which  must  have  the  opposite  sign. 
For  non-reflective  (open-end)  boundary  condition,  set  the  ghost  cell  state  identical  to  the 
boundary  cell  state.  Then  for  given  Si  at  some  instant  t,  calculate  the  cell  values  of  velocity 
Vj  =  Ui/pi  and  pressure  pi  from  formula  (18).  Zero  out  the  array  of  rates  Ri.  For  each 
facet  (pk  solve  the  Riemann  problem  using  the  given  gas  states  in  cells  wq  and  uji^  where 
ii  =  a{k,l)  and  12  =  a{k,2).  Knowing  the  gas  state  at  the  facet  calculate  the  fluxes  of 
mass,  momentum  and  energy  along  the  unit  normal  vector  n^,  multiply  by  the  facet  area 
Afc,  and  add  to  the  rate  Ri^  but  subtract  from  the  rate  according  to  the  choice  of  the 
normal  n^.  Finally,  divide  the  calculated  rates  Ri  by  volumes  Vi. 

This  procedure  defines  the  cell  rates  {Ri}  as  a  function  of  the  cell  states  {5*}  and 
therefore  the  evolution  of  the  gas  flow  can  be  calculated  by  the  explicit  Euler  scheme. 
There  is  not  much  reason  to  use  the  higher-order  Runge-Kutta  formulae  as  the  similarity 
solution  to  the  Riemann  problem  provides  time-independent  fluxes  at  facets.  The  size  of 
the  time  step  must  be  selected  according  to  the  Courant-Friedrichs-Lewy  (CFL)  condition 
which  requires  estimation  of  velocity  and  sound  speed  at  the  cell  interfaces.  Physically, 
this  condition  limits  the  time  step  to  the  value  small  enough  not  to  allow  Riemann  waves 
to  propagate  more  than  half  size  of  a  control  cell.  Note  that  the  described  Godunov-type 
scheme  is  first-order  accurate  in  space  and  time. 


5  Gas  dynamics  toolbox 

In  this  section  the  listings  of  two  reusable  MATLAB  files  required  to  simulate  blast 
propagation  in  a  fixed  domain  are  given.  Each  function  returns  a  structure  with  function 
handles. 

The  content  of  the  riemann.solver  .m  M-hle  shown  below  defines  a  MATLAB  function 
that  solves  the  Riemann  problem  using  the  Newton-Raphson  method.  The  acoustic  ap¬ 
proximation  is  taken  as  an  initial  guess  in  the  iterative  procedure.  The  function  evaluates 
the  combined  Hugoniot-Poisson  adiabatic  curve,  the  gas  state  at  the  interface  separating 
two  adjacent  cells,  and  the  complete  one-dimensional  state  of  gas  as  a  function  of  the 
coordinate  x.  The  latter  functionality  is  not  used  by  the  solver  but  provides  a  useful 
benchmark  solution  to  compare  with  numerical  results. 
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function  rs  =  riemann_solver(gciinma, tolerance, threshold) 

7,  rs  =  riemcuin_solver(gainina,tolercLnce, threshold) 

7. 

7«  Riemann  solver  for  polytropic  gas. 

7. 

7«  Parameters : 

7. 

7«  gamma  -  specific  heat  ratio 

7«  tolerance  -  convergence  tolerance  for  Newton’s  method  [default:  le-6] 

7«  threshold  -  iteration  threshold  for  Newton’s  method  [default:  10] 

7. 

7«  Function  handles: 

7o 

7«  rl  =  rs.post_state_density(rO,pO,pl) 

7. 

7«  Calculates  the  post-state  density  rl  given  the  pre-state  density  rO, 

7«  pre-state  pressure  pO,  and  post-state  pressure  pi. 

7. 

7«  vl  =  rs.post_state_velocity(rO,pO,vO,pl) 

7o 

7«  Calculates  the  post-state  velocity  vl  given  the  pre-state  density  rO, 

%  pre-state  pressure  pO,  pre-state  velocity  vO,  cuid  post-state  pressure 

7.  pi . 

7. 

%  [r,p,v]  =  rs.state(x,rl,pl,vl,r2,p2,v2) 

7. 

%  Calculates  density  r,  pressure  p,  velocity  v  at  coordinate  x  pointing 

%  from  state  1  to  state  2  separated  by  the  interface  x  =  0. 

7. 

%  [r ,a,p,vn,vt]  =  rs.interface_state(n,rl,al,pl,vl,r2,a2,p2,v2) 

7o 

%  Calculates  density  r,  sound  speed  a,  pressure  p,  normal  velocity  vn, 

%  tcuigential  velocity  vt)  at  the  interface  where  the  unit  normal  vector  n 

%  points  from  state  1  to  state  2. 

7o 

%  Examples : 

7o 

%  ***  Hugoniot-Poisson  adiabatic  curves  *** 

7. 

%  K  =  3;  N  =  1000;  gamma  =  linspace(l .5,2.5, K) ;  pO  =  1 ;  rO  =  1;  vO  =  0; 

7«  p  =  linspace(0,5*p0,N) ;  r  =  zeros(N,K);  v  =  zeros(N,K); 

%  info  =  cell(K,l); 

%  for  k  =  1:K 

%  rs  =  riemann_solver (gamma(k) ) ; 

%  info^k}  =  sprintf  ( ’  Wgamma  =  7og  \  gamma  (k) ) ; 

%  f  or  i  =  1 :  N 

%  r(i,k)  =  rs.post_state_density(r0,p0,p(i)) ; 

%  v(i,k)  =  rs.post_state_velocity(r0,p0,v0,p(i)) ; 

%  end 

%  end 

%  subplot (2 , 1 , 1) ;  plot(p,r);  xlabel( ’Pressure’ ) ;  ylabel(’Density’) ; 

%  legend (info, ’Location’ , ’Best’) ; 

%  subplot (2 , 1 , 2) ;  plot(p,v);  xlabel( ’Pressure ’) ;  ylabel ( ’Velocity ’) ; 

%  legend (info, ’Location’ , ’Best’) ; 

7. 

%  ***  Sod’s  shock  tube  *** 

7. 

%  rs  =  riemann_solver (1 .4) ; 

7o  [x,t]  =  ndgrid(linspace (-2 , 2 , 200)  , linspace (0 , 1 , 10) )  ; 

%  [r,p,v]  =  rs . state (x . /t , 1 , 1 , 0,0 . 125 ,0 . 1 , 0) ; 

%  figure;  waterf all (x. ’ ,t . ’ ,r . ’ ) ;  view(30,40);  title( ’Density’ ) ; 

%  figure;  waterfall (x. ’ ,t. ’ ,p. ’) ;  view(30,40);  title(’Pressure’) ; 

7o  figure;  waterf  all  (x.  ’  ,t .’ ,v. ’)  ;  view(30,40);  title( ’Velocity ’)  ; 

7. 

%  See  also  gas_dynamics_solver 

7o 

%  Copyright  2007-2008  Defence  Science  and  Technology  Organisation 
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"/X  Check  input 
if  nargin  <  3 

threshold  =  10; 
if  nargin  <  2 

tolerance  =  l.Oe-6; 

end 

end 

if  threshold  <  1 

error ( ^ Iteration  threshold  must  be  positive:  Zd’ , threshold) ; 

end 

if  tolercuice  <  l.Oe-14 

error ( ^Convergence  tolerance  must  be  positive:  '/g’ , tolerance) ; 

end 

if  gamma  <=  1.0 

error(^Specific  heat  ratio  must  be  greater  than  1:  Zg’, gamma); 

end 

7//0  Pre-calculate  constcuits 
gammaO  =  1.0/gamma; 
gammal  =  0 . 5* (gamma+1 . 0) ; 
gainma2  =  0 . 5*  (gamma-1 . 0) ; 
gammaS  =  gamma2/gammal ; 
gainma4  =  1 . 0/gamma2 ; 
gammas  =  gamma2*gamma0 ; 

/s'/o  Define  function  handles 

function  rl  =  density(r0,p0,pl) 
if  pi  >=  pO 

rl  =  rO*  (pl+gainma3*p0)  /  (gamma3*pl+p0) ; 
elseif  pi  >  0.0 

rl  =  rO*  (pl/pO) '■gammaO ; 

else 

rl  =  0.0; 

end 

end 

function  vl  =  velocity(rO,pO, vO,pl) 
if  pi  >=  pO 

vl  =  v0+(pl-p0)/sqrt(r0*(gammal*pl+gainma2*p0)) ; 
elseif  pi  >  0.0 

vl  =  v0+gamma4*sqrt  (gainma*pO/rO)  *  ( (pl/pO) '■gammaS-l .  0) ; 

else 

vl  =  v0-gamma4*sqrt (gamma*pO/rO) ; 

end 

end 

function  [vl ,vlprime]  =  velocity_ljet (rO,aO,pO,vO,pl) 
if  pi  >=  pO 

tl  =  gammal*pl+gamma2*p0 ; 
t2  =  sqrt(rO*tl); 
vl  =  v0+(pl-p0)/t2; 

vlprime  =  (1 . 0-0 . 5*gammal* (pl-pO) /tl)/t2 ; 
elseif  pi  >  0.0 

tl  =  (pl/pO) '■gammas ; 

vl  =  v0+a0*(tl-1.0)*gainma4; 

vlprime  =  gainmaO*aO*tl/pl ; 

else 

vl  =  v0-gamma4*a0; 
vlprime  =  inf; 

end 

end 

function  [r , a , p , v]  =  wave (x , rO , aO , pO , vO , r 1 , al , pi , vl ) 
if  pi  >  pO 

s  =  (rl*vl-rO*vO)/(rl-rO) ; 
if  X  >  s 
r  =  rO; 
a  =  aO; 
p  =  pO; 
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V  =  vO; 
elseif  X  <  s 

r  =  rl; 
a  =  al ; 
p  =  pi; 

V  =  vl ; 

else 

r  =  0.5*(r0+rl) ; 
p  =  0.5* (pO+pl) ; 

V  =  0.5*(v0+vl) ; 

a  =  sqrt  (gainma*p/r) ; 

end 

else 

if  X  >=  vO+aO 
r  =  rO; 
a  =  aO; 
p  =  pO; 

V  =  vO; 

elseif  X  <=  vl+al 
r  =  rl; 
a  =  al ; 
p  =  pi; 

V  =  vl ; 

else 

a  =  aO-gammaS* (vO+aO“x) ; 
r  =  rO*(a/aO)''gainma4; 
p  =  pO*  (r/rO) ''gamma; 

V  =  x-a; 

end 

end 

end 

function  [p , v]  =  newton_raphson(rl , al ,pl , vl ,r2 ,a2 ,p2 ,v2) 
wl  =  1.0/ (rl*al) ; 
w2  =  1.0/(r2*a2) ; 

p  =  max((wl*pl+w2*p2+vl-v2)/(wl+w2) .tolerance) ; 
i  =  0; 

e  =  1 . 0+tolerance ; 
while  e  >  tolerance 
i  =  i+1; 

if  i  >  threshold 

error ( 'Reached  iteration  threshold:  /og’.e); 

end 

[ujuprime]  =  velocity_ljet(rl,al,pl,-vl,p) ; 

[v.vprime]  =  velocity_ljet(r2,a2,p2,v2,p) ; 
e  =  (v+u)/(vprime+uprime) ; 
p  =  max (p-e, tolerance) ; 
e  =  abs(e/p); 

end 

V  =  0.5*(v-u) ; 

end 

function  [r,p,v]  =  State(x,rl,pl,vl,r2,p2,v2) 
al  =  sqrt  (gainma*pl/rl) ; 
a2  =  sqrt  (gainma*p2/r2) ; 

[pO,vO]  =  newton_raphson(rl,al,pl,vl,r2,a2,p2,v2) ; 
r  =  zeros (size (x) ) ; 
p  =  zeros (size (x) ) ; 

V  =  zeros (size (x) ) ; 
for  i  =  l:numel(x) 

if  x(i)  <  vO 

rO  =  density(rl ,pl ,p0) ; 
aO  =  sqrt  (gainma*pO/rO) ; 

[r(i) ,a,p(i) ,v(i)]  =  wave (-x(i) ,rl ,al ,pl ,-vl ,r0, aO,pO, -vO) ; 
v(i)  =  -v(i); 
elseif  x(i)  >  vO 

rO  =  density(r2 ,p2 ,pO) ; 
aO  =  sqrt  (gainma*pO/rO) ; 

[r(i) ,a,p(i) ,v(i)]  =  wave(x(i) ,r2,a2,p2,v2,r0,a0,p0,v0) ; 
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else 

rO  =  densityCrl ,pl ,pO) ; 
aO  =  sqrt  (gainma*pO/rO) ; 

[rrl,a,ppl,vvl]  =  wave(-x(i) ,rl,al,pl,-vl,rO,aO,pO,-vO) ; 
vvl  =  -vvl; 

rO  =  density(r2,p2,p0) ; 
aO  =  sqrt  (gainma*pO/rO) ; 

[rr2,a,pp2, vv2]  =  wave(x(i) ,r2,a2,p2,v2,r0,a0,p0,v0) ; 
r(i)  =  0.5*(rrl+rr2) ; 
p(i)  =  0.5*(ppl+pp2) ; 
v(i)  =  0.5*(vvl+vv2); 

end 

end 

end 

function  [r , a , p , vn , vt]  =  interf ace_state (n , r 1 , al ,pl , vl , r2 , a2 ,p2 , v2) 
vln  =  dot (vl ,n) ; 
v2n  =  dot (v2 ,n) ; 
vlt  =  vl-vln*n; 
v2t  =  v2-v2n*n; 

[pO,vO]  =  newton_raphson(rl,al,pl,vln,r2,a2,p2,v2n) ; 
switch  sign(vO) 
case  1 

rO  =  density(rl ,pl ,p0) ; 
aO  =  sqrt  (gainma*pO/rO) ; 

[r,a,p,vn]  =  wave(0.0,rl,al,pl,-vln,r0,a0,p0,“v0) ; 
vn  =  -vn; 
vt  =  vlt; 
case  -1 

rO  =  density(r2 ,p2 ,p0) ; 
aO  =  sqrt  (gainma*pO/rO) ; 

[r,a,p,vn]  =  wave (0 . 0,r2, a2 ,p2 , v2n,r0 ,a0 ,p0 , vO) ; 
vt  =  v2t; 
otherwise 

rO  =  density(rl ,pl ,p0) ; 
aO  =  sqrt  (gainma*pO/rO) ; 

[rrl , aal ,ppl ,vnl]  =  wave(0.0,rl,al,pl,-vln,r0,a0,p0,-v0) ; 
vnl  =  -vnl ; 

rO  =  density(r2,p2,p0) ; 
aO  =  sqrt  (gainma*pO/rO) ; 

[rr2,aa2,pp2,vn2]  =  wave(0.0,r2,a2,p2,v2n,r0,a0,p0,v0) ; 
r  =  0.5*(rrl+rr2) ; 
p  =  0 . 5* (ppl+pp2) ; 
vn  =  0 . 5* (vnl+vn2) ; 
vt  =  0 . 5* (vlt+v2t) ; 
a  =  sqrt  (gaiiima*p/r) ; 

end 

end 

"/X  Publish  the  interface 
rs. state  =  Estate; 

rs . interf ace_state  =  Ointerf ace_state; 
rs .post_state_density  =  ©density; 
rs .post_state_velocity  =  ©velocity; 
end 


The  content  of  the  gas_dynamics_solver  .m  M-file  shown  below  defines  a  MATLAB 
function  that  implements  a  gas  dynamics  solver.  The  user  needs  to  populate  a  data 
structure  and  pass  it  to  the  solver.  The  solver  is  generic  as  it  does  not  use  any  specific 
mesh.  Instead,  only  a  facet-to-cell  connectivity  matrix,  an  array  of  facet  vector  area 
elements  and  an  array  of  control  volumes  are  used  to  define  geometry.  Moreover,  this 
function  does  not  depend  on  the  riemann.solver  .m  function  as  it  works  with  a  function 
handle  for  interface  state  calculation.  The  riemann_solver  .m  provides  such  a  handle,  but 
other  implementations  can  be  equally  used,  such  as  the  approximate  Roe  solver  [Roe  1981]. 
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function  gds  =  gas_dynainics_solver 
7,  gds  =  gas_dynainics_solver 

7. 

7«  Gas  dynamics  solver. 

7. 

7«  Function  handles: 

7. 

7«  model  =  gds .  initialise  (model) 

7o 

7«  Validates  the  model  structure,  allocates  auxiliary  fields  and  populates 
7«  the  conserved  variables . 

7o 

7«  model  =  gds. update_rates (model) 

7. 

7«  Calculates  the  interface  states  from  the  current  physical  variables  cuid 
7«  updates  the  rates  of  the  conserved  variables. 

7o 

7«  model  =  gds. update_states (model) 

7. 

7«  Updates  the  physical  variables  from  the  conserved  variables. 

7«  Should  be  always  called  after  advancing  the  time  step  cuid  appying 
%  boundary  conditions  that  must  affect  the  conserved  variables  only. 

7. 

%  The  model  structure  must  contain  the  following  fields: 

7o 

7«  model. gamma  -  polytropic  index 

%  model. state  -  interface  state  sampling  handle 

%  model. fc  -  f acet-to-cell  connectivity  matrix 

%  model. ds  -  vector  area  elements 

7  model. dv  -  volume  elements 

%  model. r  -  density 

%  model. p  -  pressure 

7«  model. V  -  velocity 

7o 

%  The  following  auxiliary  fields  will  be  allocated: 

7. 

%  model. gammal  -  model. gamma  -  1 

%  model. da  -  facet  area  elements 

%  model. n  -  facet  unit  normal  vectors 

%  model. a  -  sound  speed 

%  model. e  -  total  energy 

7  model. u  -  momentum 

7  model. r_rate  -  density  rate 

7  model. e_rate  -  total  energy  rate 
7  model. u_rate  -  momentum  rate 

7  model. ir  -  interface  density 

7  model. ip  -  interface  pressure 

7  model. ia  -  interface  sound  speed 

7  model. ivn  -  interface  normal  component  of  velocity 

7 

7  Examples : 

7 

7  ***  Sod's  shock  tube  *** 

7 

7  model. gamma  =  1.4;  rs  =  r iemann_solver (model. gamma) ; 

7  model. state  =  rs . interf ace_state ;  gds  =  gas_dynamics_solver ; 

7  N  =  200;  L  =  1.0;  dx  =  L/N;  x  =  -dx/2 : dx : l+dx/2; 

7  model. dv  =  dx*ones (N+2 , 1) ;  model. ds  =  ones(N+l,l); 

7  model. fc  =  [1 :N+1 ; 2 :N+2] . ' ;  model. r  =  zeros (N+2 , 1) ; 

7  model. p  =  zeros (N+2 , 1) ;  model. v  =  zeros (N+2 , 1) ; 

7  boundary_cells  =  [2,N+1];  ghost_cells  =  [l,N+2]; 

7  for  i  =  l:N+2 

7  if  x(i)  <  L/2 

7  model. r(i)  =  1.0;  model. p(i)  =  1.0; 

7  else 

7  model. r(i)  =  0.125;  model. p(i)  =  0.1; 

7  end 


13 


DSTO-TR-2139 


7,  end 

7«  model  =  gds . initialise(model) ; 

7o  t  =  0;  tf  =  0.25; 

7o  while  t  <  tf 

7«  model  =  gds. update_rates (model) ; 

7«  dt  =  dx/max  (model .  ia+abs  (model .  ivn) ) ; 

7o  dt  =  min(dt,tf-t)  ;  t  =  t+dt; 

7«  model,  r  =  model.  r+model.r_rate*dt; 

7«  model,  e  =  model .  e+model .  e_rate*dt ; 

7«  model,  u  =  model.  u+model.u_rate*dt; 

7«  model. r(ghost_cells)  =  model. r (boundary_cells) ; 

7«  model .  e  (ghost_cells)  =  model .  e  (boundary_cells) ; 

7«  model  .u(ghost_cells)  =  model  .u(boundary_cells) ; 

7«  model  =  gds. update_states (model) ; 

7o  subplot  (3 , 1 , 1)  ; 

7o  plot  (x(2  :  end-1)  .model  .r  (2  :  end-1)  ,  ^  .  ’)  ; 

7«  title (sprintf ( 'Density  at  t  =  7og’,t)); 

7o  subplot  (3 , 1 , 2)  ; 

7o  plot  (x(2  :  end-1)  .model  .p (2  :  end-1)  .'.'); 

7«  title (sprintf ( 'Pressure  at  t  =  7g',t)); 

7o  subplot (3 , 1 ,3)  ; 

%  plot (x(2 : end-1) .model .v(2 : end-1) .'.'); 

%  title  (sprintf  ( 'Velocity  at  t  =  ^gSt)); 

7  drawnow; 

%  end 

% 

%  See  also  riemann_solver 

% 

%  Copyright  2007-2008  Defence  Science  and  Technology  Organisation 

VL  Define  function  handles 

function  model  =  initialise(model) 

7«  Check  required  fields 
if  ''isfield(model, 'gamma') 

error ( 'Undefined  polytropic  index:  model . gamma' ) ; 

end 

if  ''isfield(model .’ state  ’ ) 

error ( 'Undefined  interface  state  sampling  handle:  model . state ’) ; 

end 

if  "isa (model .  state .  ’  function_hcuidle ' ) 

error(’Field  model. state  must  be  a  function  handle’); 

end 

if  "isf ield(model . ’ f c ’ ) 

error ( 'Undefined  f acet-to-cell  connectivity  matrix:  model. fc'); 

end 

if  ''isfield(model. 'ds') 

error ( 'Undefined  vector  area  elements:  model. ds’); 

end 

if  ''isfield(model. 'dv') 

error ( 'Undefined  volume  elements:  model. dv’); 

end 

if  ''isfield(model. 'r') 

error ( 'Undefined  density:  model. r’); 

end 

if  ''isfield(model. 'p') 

error ( 'Undefined  pressure:  model. p’); 

end 

if  ''isfield(model, 'v') 

error ( 'Undefined  velocity:  model. v’); 

end 

7#  Check  consistency 

if  size (model . ds . 2)  ~=  size (model . v. 2) 

error ( 'Unequal  number  of  columns  of  model. ds  and  model. v’); 

end 

model. DIM  =  size(model.v,2) ; 
if  model. DIM  <  1  I  I  model. DIM  >  3 

error(  ’  Invalid  space  dimension:  7od’ .model  .DIM) ; 
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end 

if  size (model . fc , 2)  ~=  2 

error( ’Connectivity  matrix  must  have  two  columns’); 

end 

if  size(model.dv,2)  ~=  1 

error( ’Control  volumes  must  form  a  column  vector’); 

end 

if  size(model.r,2)  "=  1 

error( ’Density  field  must  be  a  column  vector’); 

end 

if  size(model.p,2)  "=  1 

error( ’Pressure  field  must  be  a  column  vector’); 

end 

if  size (model . dv, 1)  ~=  size(model.r,l) 

error ( ’Unequal  number  of  rows  of  model. dv  and  model.r’); 

end 

if  size (model . dv, 1)  ~=  size(model.p,l) 

error ( ’Unequal  number  of  rows  of  model. dv  and  model.p’); 

end 

if  size(model.dv,l)  ~=  size (model . v, 1) 

error ( ’Unequal  number  of  rows  of  model. dv  and  model. v’); 

end 

7o  Initialise 

model. gammal  =  model.gamma*l .0; 
model. da  =  sqrt (dot (model. ds, model. ds, 2) ) ; 
model. n(:,l)  =  model . ds (:, 1) . /model .da; 
if  model. DIM  >  1 

model. n(:, 2)  =  model .ds (: ,2) . /model . da; 
if  model. DIM  >  2 

model. n(:, 3)  =  model.ds(: ,3) ./model.da; 

end 

end 

model. a  =  zeros(size(model.r)) ; 
model. e  =  zeros(size(model.p)) ; 
model. u  =  zeros(size(model.v)) ; 
model. r_rate  =  zeros(size(model.r)) ; 
model. e_rate  =  zeros(size(model.e)) ; 
model. u_rate  =  zeros(size(model.u)) ; 
model. ir  =  zeros (size (model . ds , 1) , 1) ; 
model. ip  =  zeros (size (model . ds , 1) , 1) ; 
model. ia  =  zeros (size (model . ds , 1) , 1) ; 
model. ivn  =  zeros (size (model . ds , 1) , 1) ; 
model,  a  =  sqrt  (model.  gainma*model.p. /model,  r) ; 
model. u(:,l)  =  model. r. *model. v(: ,1) ; 
if  model. DIM  >  1 

model. u(:, 2)  =  model .r . *model . v( : ,2) ; 
if  model. DIM  >  2 

model. u(:, 3)  =  model.r.*model.v(: ,3) ; 

end 

end 

model. e  =  model. p/model. gainmal+0.5*dot(model.u, model. v, 2) ; 

end 

function  model  =  update_rates (model) 

model. r_rate  =  zeros(size(model.r_rate)) ; 
model. e_rate  =  zeros(size(model.e_rate)) ; 
model. u_rate  =  zeros(size(model.u_rate)) ; 
for  k  =  1 : size (model . fc , 1) 

11  =  model. fc(k, 1) ; 

12  =  model.fc(k,2) ; 

[r,a,p,vn,vt]  =  model. state(model. n(k, :),.. . 

model . r (il) , model . a(il) ,model .p(il) ,model . v(il ,:),... 
model . r (i2) , model . a(i2) ,model .p(i2) ,model . v(i2 , : ) ) ; 
model. ir(k)  =  r; 
model. ia(k)  =  a; 
model. ip(k)  =  p; 
model. ivn (k)  =  vn; 

V  =  vt+vn*model .n(k, : ) ; 
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r_flux  =  r*vn*model . da(k) ; 

e_flux  =  r_flux*(a*a/inodel.gainmal+0.5*dot(v,v,2)) ; 
u_flux  =  r_flux*v+p*model.ds(k, :) ; 
model . r_rate (il)  =  model . r_rate (il) -r_f lux; 
model . r_rate (i2)  =  model . r_rate (i2)+r_f lux; 
model . e_rate (il)  =  model . e_rate (il) -e_f lux; 
model . e_rate (i2)  =  model . e_rate (i2)+e_f lux; 
model. u_rate ( i 1, :)  =  model .u_rate(il )-u_f lux; 
model. u_rate(i2, :)  =  model .u_rate(i2 )+u_f lux; 

end 

model. r_rate  =  model. r_rate. /model. dv; 
model. e_rate  =  model. e_rate. /model. dv; 
model .u_rate( 1)  =  model .u_rate( 1) . /model . dv; 
if  model. DIM  >  1 

model.u_rate( : ,2)  =  model.u_rate(: ,2) ./model.dv; 
if  model. DIM  >  2 

model .u_rate( : ,3)  =  model .u_rate( : ,3) . /model .dv; 

end 

end 

end 

function  model  =  update_states (model) 

model. v(:,l)  =  model. u(: ,1) ./model. r; 
if  model. DIM  >  1 

model. v(:, 2)  =  model .u( : ,2) . /model .r ; 
if  model. DIM  >  2 

model. v(:, 3)  =  model. u(: ,3) ./model. r; 

end 

end 

model. p  =  model. gammal*(model.e-0.5*dot(model.u, model. V, 2)) ; 
model. a  =  sqrt (model. gamma*model.p. /model. r) ; 

end 

7//o  Publish  the  interface 
gds . initialise  =  ©initialise; 
gds.update_rates  =  ©update_rates ; 
gds.update_states  =  ©update_states ; 
end 


6  Usage  examples 


6.1  Sod’s  shock  tube 


The  exact  similarity  solution  to  Sod’s  shock  tube  [Sod  1978]  specified  by  the  following 
initial  piecewise-constant  state  of  gas 


(1.0,  1.0,  0.0)  {x  <  0) 
(0.125,  0.1,  0.0)  (x  >  0) 


(40) 


provides  a  useful  benchmark  for  gas  dynamics  solvers. 

The  following  MATLAB  script  simulates  the  benchmark  problem  of  ideal  polytropic 
gas  with  specific  heat  ratio  7  =  1.4  on  a  uniform  mesh  consisting  of  400  control  cells.  The 
numerical  results  are  displayed  in  Figure  2.  It  is  seen  that  the  numerical  solution  and  the 
exact  benchmark  solution  are  in  good  agreement.  Note  that,  due  to  numerical  viscosity, 
the  contact  discontinuity  is  smeared  out  to  some  extent  whereas  the  shock  front  position 
is  captured  quite  accurately. 
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(i)  Pressure.  (ii)  Density. 


(iii)  Velocity. 


(iv)  Internal  energy. 


Figure  2:  Gas  state  in  Sod’s  shock  tube  at  time  t  =  0.25. 


Xk  Sod’s  shock  tube  simulation 
xl  =  -0.5;  rl  =  1.0;  pi  =  1.0;  vl  =  0.0; 
x2  =  0.5;  r2  =  0.125;  p2  =  0.1;  v2  =  0.0; 
tf  =  0.25;  NX  =  400; 

Xk  Benchmark  solution 

model. gamma  =  1.4;  rs  =  riemann_solver (model. gamma) ; 
model. state  =  rs . interf ace_state ; 

X  =  linspace(xl ,x2 , 1000) ;  [R,P,V]  =  rs . state (X/tf ,rl , pi ,vl ,r2 ,p2 ,v2) ; 
Xk  Generate  a  uniform  mesh 

dx  =  (x2-xl)/NX;  x  =  (xl-dx/2:dx:x2+dx/2) . ’ ; 

model. fc  =  [1 : NX+1 ; 2 :NX+2]  . ’ ;  model. ds  =  ones (NX+1 , 1) ; 

model. dv  =  dx*ones (NX+2 , 1) ;  B  =  [2, NX+1];  G  =  [l,NX+2];  I  =  2:NX; 

Xk  Run  simulation 

disp(’***  Sod’’s  shock  tube  ♦**’); 

t  =  0.0; 

[model. r, model. p, model. v]  =  rs.state(x/t,rl,pl,vl,r2,p2,v2) ; 
model. V  =  zeros (NX+2 , 1) ;  gds  =  gas_dynamics_solver ; 
model  =  gds . initialise (model) ; 
while  t  <  tf 

model  =  gds .update_rates(model) ; 
dt  =  dx/max (model . ia+abs (model . ivn) ) ; 
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dt  =  minCdt  ,tf-t)  ;  t  =  t+dt;  dispCsprintf  ( ’t  =  Vog^jt)); 

model. r  =  model . r+model . r_rate*dt ; 

model. e  =  model . e+model . e_rate*dt ; 

model. u  =  model. u+model.u_rate*dt; 

model. r(G)  =  model. r(B); 

model. e(G)  =  model. e(B); 

model. u(G)  =  -model. u(B); 

model  =  gds .update_states (model) ; 

end 


7//0  Plot  the  final  state 

plot(X,R» ’r-’ ,x, model. r, ^b. O ;  grid  on; 

axis ( [xl ,x2 ,0, 1 . 2] ) ;  legend ( ’ exact ’ , 'numerical ’ , ’Location’ , ’Best ’ ) ; 
print  “dpng  sod_shock_tube_density ; 

plot(X,P./R/model.gainmal, ’r-’ ,x, model. p. /model. r /model. gammal, ’b. ’) ;  grid  on; 
axis ( [xl ,x2 , 1 . 0,3 . 0] ) ;  legend ( ’ exact ’ , ’numerical ’ , ’Location’ , ’Best ’ ) ; 
print  -dpng  sod_shock_tube_energy ; 
plot(X,P, ’r-’ ,x, model. p, ’b. ’) ;  grid  on; 

axis ( [xl ,x2 ,0, 1 . 2] ) ;  legend ( ’ exact ’ , ’numerical ’ , ’Location’ , ’Best ’ ) ; 
print  -dpng  sod_shock_tube_pressure ; 
plot(X,V, ’r-’ ,x, model. V, ’b. ’) ;  grid  on; 

axis ( [xl ,x2 ,0, 1 . 2] ) ;  legend ( ’ exact ’ , ’numerical ’ , ’Location’ , ’Best ’ ) ; 
print  -dpng  sod_shock_tube_velocity; 


6.2  Rotationally  symmetric  blast 


Another  interesting  one-dimensional  example  that  can  be  used  as  a  benchmark  for 
testing  multi-dimensional  code  is  a  rotationally  symmetric  blast  in  two  or  three  dimensions. 
All  the  gas  state  variables  are  functions  of  the  radial  coordinate  r  and  time  t  and  velocity 
has  only  the  radial  non-zero  component  v.  The  governing  equations  are  as  follows 


^  d{pv)  ^  ^ 

dt  dr  r  ’ 


(41) 


d{pv)  d[[pv‘^+p)v\  V  pv"^ 

dt  dr  r 


ds  9[(e  +p)  v]  ly  {e  +  p)  V 
dt  dr  r 


(42) 

(43) 


where  u  =  2  for  rotationally  symmetric  three-dimensional  blast  (described  in  the  spherical 
coordinates)  and  v  =  1  for  rotationally  symmetric  two-dimensional  blast  (described  in  the 
polar  coordinates).  Note  that  the  case  ly  =  0  corresponds  to  the  one-dimensional  problem 
describing  shock  propagation  in  a  tube,  such  as  the  one  considered  in  the  previous  section. 
For  nonzero  ly,  equations  (41)-(43)  contain  geometric  source  terms  which  can  be  easily 
incorporated  into  the  numerical  scheme  after  the  calculation  of  the  gas  state  rates  due  to 
fluxes  alone. 


The  following  script  simulates  the  propagation  of  blast  waves  from  a  spherical  charge 
(ly  =  2).  The  results  of  the  numerical  simulation  are  presented  in  Figure  3.  It  is  seen  the 
formation  of  a  secondary  shock  wave  following  the  main  shock  [Erode  1959]. 


/s'/o  Spherically  symmetric  blast 

clear;  project_name  =  ’radial_explosion’ ;  choice  =  1; 
if  exist  (sprintf  (’ /oS  .mat  ’  ,project_name) ,  ’file’) 

choice  =  menu ( ’Re-run? ’Yes ’No,  just  post-process’ , ’Cancel’) ; 
pause(0 . 1) ; 
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(i)  Pressure.  (ii)  Density. 


(iii)  Velocity. 


(iv)  Temperature. 


Figure  3:  Blast  waves  from  a  spherical  charge. 


end 

switch  choice 
case  1 

7,  Specify  parameters 

model. gamma  =  1.4;  rs  =  riemcuin_solver (model . gamma, 1 . Oe-9 , 100) ; 

model. state  =  rs . interf ace_state ; 

nu  =  2;  L  =  100.0;  w  =  1.0;  R  =  287.0; 

pO  =  1.01e5;  rO  =  1.204;  pi  =  1.0e5*p0;  rl  =  l.OeS; 

NX  =  1000;  CFL  =0.9; 

7o  Generate  a  uniform  mesh 

dx  =  L/NX;  X  =  (-dx/2:dx:L+dx/2) . ^ ;  z  =  -nu./x; 
model. dv  =  dx*ones (NX+2 , 1) ;  model. ds  =  ones (NX+1 , 1) ; 
model. fc  =  [1 : NX+1 ; 2 :NX+2] . ' ; 

B1  =  2;  G1  =  1;  B2  =  NX+1;  G2  =  NX+2; 

7«  Initial  conditions 
model. r  =  zeros (NX+2 , 1) ; 
model. p  =  zeros (NX+2 , 1) ; 
model. V  =  zeros (NX+2 , 1) ; 
for  i  =  l:NX+2 
if  x(i)  <  w 

model. r(i)  =  rl;  model. p(i)  =  pi; 

else 

model. r(i)  =  rO;  model. p(i)  =  pO; 

end 
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end 

7o  Allocate  variables  for  post-processing 

times  =  [0.025,0.05,0.075,0.1,0.125,0.15];  K  =  length (times) ; 

pressures  =  zeros(NX,K); 

velocities  =  zeros(NX,K); 

densities  =  zeros(NX,K); 

temperatures  =  zeros(NX,K); 

I  =  2:NX+1;  x  =  x(I) ; 

7o  Run  simulation 

disp(^***  Spherically  symmetric  blast 

gs  =  gas_dynamics_solver ;  model  =  gs . initialise (model) ; 
t  =  0.0; 
for  k  =  1:K 

tf  =  times (k) ; 
while  t  <  tf 

model  =  gs. update_rates (model) ; 
dt  =  CFL*dx/max(abs (model . ivn)+model . ia) ; 
dt  =  min(dt,tf-t) ;  t  =  t+dt; 
model. r  =  model . r+model . r_rate*dt ; 
model. e  =  model . e+model . e_rate*dt ; 
model. u  =  model. u+model.u_rate*dt; 
model  =  gs. update_states (model) ; 
model. r_rate  =  z.*model.u; 

model. e_rate  =  model. r_rate.* (model. e+model.p) ./model. r 

model. u_rate  =  model. r_rate.*model.v; 

model. r  =  model . r+model . r_rate*dt ; 

model. e  =  model . e+model . e_rate*dt ; 

model. u  =  model. u+model.u_rate*dt; 

model. r(Gl)  =  model. r(Bl);  %  reflective  BC 

model. e(Gl)  =  model. e(Bl); 

model. u(Gl)  =  -model .u(Bl) ; 

model. r(G2)  =  model. r(B2);  %  open-end  BC 

model. e(G2)  =  model. e(B2); 

model. u(G2)  =  model. u(B2); 

model  =  gs. update_states (model) ; 

plot(x,model.p(I) , ’k. ’) ; 

title(sprintf ('Pressure  at  t  =  7og',t)); 

drawnow ; 

end 

pressures(: ,k)  =  model. p(I); 
velocities(: ,k)  =  model. v(I); 
densities( : ,k)  =  model. r(I); 

temperatures (: ,k)  =  model. p(I) ./model. r(I)/R; 

end 

save(project_name) ; 
case  2 


load(pro ject_name) ; 
otherwise 
return; 

end 

7«7o  Generate  images 
info  =  cell (size (times) ) ; 
for  k  =  1:K 

info-Ck}  =  sprintf('t  =  7»g  ms  '  ,  1000 . 0*times  (k)  )  ; 

end 

figure;  plot(x,pressures/p0) ;  grid  on; 

xlabeK 'Distance  (m)’);  ylabeK’ Pressure  (atm)'); 

legend (info, 'Location' , 'Best') ; 

eval (sprintf  ( 'print  -dpng  7oS_pressure’ ,project_name)) ; 
figure;  plot(x,velocities) ;  grid  on; 
xlabeK 'Distance  (m)’);  ylabeK’Velocity  (m/s)’); 
legend (info, 'Location' , 'Best') ; 

eval  (sprintf  ( 'print  -dpng  7oS_velocity  ’  ,project_name) ) ; 
figure;  plot(x,densities) ;  grid  on; 
xlabeK 'Distance  (m)');  ylabeK 'Density  (kg/m''3)’); 
legend (info, 'Location' , 'Best') ; 

evaKsprintf  ( 'print  -dpng  7oS_density' ,project_name)) ; 
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figure;  plot(x,teinperatures) ;  grid  on; 

xlabeK ^Distance  (m)’);  ylabeK’ Temperature  (K)  O  ; 

legend (info, ’Location^ , ’Best’) ; 

eval (sprintf  ( ’print  -dpng  /(s.temperature’ ,project_naine)) ; 


6.3  Two-dimensional  blast  propagation 

Consider  a  two-dimensional  explosion  in  air  occupying  a  square  domain  with  a  side 
length  of  10  m  at  initial  atmospheric  pressure  of  101  kPa  (1  atmosphere)  and  temperature 
of  20°C.  The  blast  is  modelled  by  a  localised  region  of  10-atmosphere  over-pressure  and 
temperature  of  2000°C. 

The  following  MATLAB  script  simulates  this  problem  on  a  uniform  mesh.  The  pressure 
and  temperature  fields  for  increasing  times  are  presented  in  Figures  4  and  5  respectively. 


"/X  2D  blast  propagation  in  a  square  domain 
project_name  =  ’blast2d’; 

disp (sprintf  ( ’Project  name:  7,s’ ,project_name)) ; 
gds  =  gas_dynamics_solver ; 

R  =  287.0;  KelvinO  =  273.0; 

TO  =  KelvinO+20.0;  pO  =  l.OleS;  rO  =  pO/TO/R; 

T1  =  KelvinO+2000.0;  pi  =  ll*p0;  rl  =  pl/Tl/R; 

CFL  =  0.5;  save_time  =  0.002;  kmax  =  6; 

7//0  Pre-process  if  necessary 
k  =  0; 

if  "exist  (sprintf  ( ’7os7od.mat  ’  ,project_name ,k) ,  ’ f  ile ’ ) 

7o  Generate  mesh 

s.LX  =  10.0;  s.LY  =  10.0;  s.NX  =  250;  s.NY  =  250; 
s.map  =  @(ix,iy)  sub2ind( [s .NX+2 , s . NY+2] , ix, iy) ; 
dx  =  s.LX/s.NX;  dy  =  s.LY/s.NY;  s.dx  =  min( [dx,dy] ) ; 
s.x  =  -dx/2:dx:s.LX+dx/2;  s.y  =  -dy/2 : dy : s .LY+dy/2; 

NC  =  (s.NX+2)*(s.NY+2) ;  NF  =  (s . NX+1) *s .NY+s .NX* (s .NY+1) ; 
disp  (sprintf  (’7od  facets,  7od  cells  ’  ,NF,NC-4)  )  ; 

s.dv  =  (dx*dy) *ones (NC , 1) ;  s.fc  =  zeros(NF,2);  s.ds  =  zeros(NF,2); 

s.BE  =  zeros(s.NY,l)  ;  s.GE  =  zeros  (s  .NY,  1)  ;  7o  east  BC 

s.BW  =  zeros(s.NY,l)  ;  s.GW  =  zeros  (s  .NY,  1)  ;  7o  west  BC 

s.BS  =  zeros(s.NX,l)  ;  s.GS  =  zeros  (s  .NX ,  1)  ;  7o  south  BC 

s.BN  =  zeros(s.NX,l)  ;  s.GN  =  zeros  (s  .NX ,  1)  ;  7o  north  BC 

i  =  0;  7o  facet  index 

for  ix  =  l:s.NX+l 

for  iy  =  2:s.NY+l 

i  =  i+1;  s.ds(i,:)  =  [dy,0] ; 
s.fc(i,l)  =  s .map(ix, iy) ; 
s.fc(i,2)  =  s.map(ix+l,iy) ; 

end 

end 

for  iy  =  l:s.NY+l 

for  ix  =  2:s.NX+l 

i  =  i+1;  s.ds(i,:)  =  [0,dx] ; 
s.fc(i,l)  =  s .map(ix, iy) ; 
s.fc(i,2)  =  s .map(ix, iy+1) ; 

end 

end 

for  iy  =  2:s.NY+l 

s.BE(iy-l)  =  s.map(2,iy);  s.GE(iy-l)  =  s.map(l,iy); 
s.BW(iy-l)  =  s.map(s.NX+l,iy) ;  s.GW(iy-l)  =  s .map(s .NX+2 , iy) ; 

end 

for  ix  =  2:s.NX+l 

s.BS(ix-l)  =  s.map(ix,2);  s.GS(ix-l)  =  s.map(ix,l); 
s.BN(ix-l)  =  s.map(ix,s.NY+l) ;  s.GN(ix-l)  =  s .map (ix, s .NY+2) ; 

end 
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%  Initial  conditions 

xa=7.0;  ya=8.0;  a  =0.5;  b=  100.0; 
density  =  @(x,y) . . . 

0 . 5*  ( (rl+rO)  +  (rl-rO)  *tainh(b*  (a-sqrt  ( (x-xa)  “2+  (y-ya)  ''2) ) ) ) ; 
pressure  =  @(x,y)... 

0 . 5*  ( (pl+pO)  +  (pl-pO)  *tanh(b*  (a-sqrt  ( (x-xa)  “2+  (y-ya)  ''2) ) ) ) ; 
s.r  =  zeros(NC,l);  s.p  =  zeros(NC,l);  s.v  =  zeros(NC,2); 
for  ix  =  l:s.NX+2 

for  iy  =  l:s.NY+2 

i  =  s .map(ix, iy) ; 

s.r(i)  =  density(s.x(ix) ,s.y(iy)) ; 
s.p(i)  =  pressure(s.x(ix) ,s.y(iy)) ; 

end 

end 

%  Specify  Riemann  solver 

s. gamma  =  1.4;  rs  =  riemann_solver(s .gamma) ; 

s. state  =  rs . interf ace_state ;  s.t  =  0.0;  s  =  gds . initialise (s) 
eval (sprintf  ( ’ save  7«s%d  -struct  s’ ,project_name,k)) ; 

end 


7, 7  Find  the  latest  restart  file 

while  exist(sprintf  (’7»s7od.mat’  ,project_name ,k) ,  ’file’) 
k  =  k+1; 

end 

k  =  k-1; 

7«7o  Load  start/restart  file 
s  =  load(sprintf  (’7os7«d’ ,project_name,k)) ; 
while  k  <  kmax 

elapsed_time  =  cputime; 
k  =  k+1; 

tf  =  s.t  +  save_time; 
while  s.t  <  tf 

s  =  gds .update_rates(s) ; 

dt  =  CFL*s.dx/max(s.ia+abs(s.ivn)) ; 

dt  =  min(dt,tf-s.t) ;  s.t  =  s.t+dt; 

disp(sprintf  ( ’t  =  7og’,s.t)); 

s.r  =  s.r+s.r_rate*dt; 

s.e  =  s . e+s . e_rate*dt ; 

s.u  =  s.u+s.u_rate*dt; 

7o  reflective  east  BC 

s.r(s.GE)  =  s.r(s.BE);  s.e(s.GE)  =  s.e(s.BE); 
s.u(s.GE,l)  =  -s.u(s.BE,l) ;  s.u(s.GE,2)  =  s.u(s.BE,2); 
7o  reflective  west  BC 

s.r(s.GW)  =  s.r(s.BW);  s.e(s.GW)  =  s.e(s.BW); 
s.u(s.GW,l)  =  -s.u(s.BW,l) ;  s.u(s.GW,2)  =  s.u(s.BW,2); 
7«  reflective  south  BC 

s.r(s.GS)  =  s.r(s.BS);  s.e(s.GS)  =  s.e(s.BS); 
s.u(s.GS,l)  =  s.u(s.BS,l);  s.u(s.GS,2)  =  -s .u(s . BS , 2) ; 
7»  reflective  north  BC 

s.r(s.GN)  =  s.r(s.BN);  s.e(s.GN)  =  s.e(s.BN); 
s.u(s.GN,l)  =  s.u(s.BN,l);  s.u(s.GN,2)  =  -s .u(s . BN, 2) ; 
7«  synchronise  physical  variables 
s  =  gds .update_states (s) ; 

end 

elapsed_time  =  cputime-elapsed_time ; 

disp(sprintf (’CPU  time:  7«g  secs  /  7»g  mins  /  7«g  hours’,... 

elapsed_time , elapsed_time/60 ,elapsed_time/3600) ) ; 
eval (sprintf  ( ’ save  7«s7od  -struct  s’ ,project_name,k)) ; 

end 

7«7o  Post  process 

disp( ’Generating  PNG  images...’); 
set(clf , ’Renderer’ , ’ZBuffer’) ; 
k  =  0; 

while  exist(sprintf  (’7os7od.mat’  ,project_name,k) ,  ’file’) 
s  =  load(sprintf  (’7«s7»d’ ,project_name,k)) ; 
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[x,y]  =  ndgrid(s.x,s.y) ;  z  =  zeros (size (x) ) ; 
for  ix  =  2:s.NX+l 

for  iy  =  2:s.NY+l 

i  =  s .map(ix, iy) ; 
z(ix,iy)  =  s.p(i)/1000; 

end 

end 

surf (x(2 : end-1 , 2 : end-2) ,y(2 : end-1 , 2 : end-2) ,z(2 : end-1,2: end-2) ) ; 
title(sprintf  ( ’Pressure  (kPa)  at  time  t  =  '/og  ms  ’ ,  1000*s  .t) ) ; 
axis (  [0 , s . LX,0 , s .LY, 0 ,pl/2000] ) ;  shading  interp; 
evaKsprintf  (’print  -dpng  ‘/.s'/od.pressure’ ,project_name,k)) ; 
for  ix  =  2:s.NX+l 

for  iy  =  2:s.NY+l 

i  =  s .map(ix, iy) ; 
z(ix,iy)  =  s.p(i)/s.r(i)/R; 

end 

end 

surf (x(2 : end-1 , 2 : end-2) ,y(2 : end-1 , 2 : end-2) ,z(2 : end-1,2: end-2) ) ; 
title(sprintf  ( ’Temperature  (K)  at  time  t  =  7og  ms’ ,1000*s.t)) ; 
axis ( [0 , s . LX,0, s .LY, 0 ,T1] ) ;  shading  interp; 

evaKsprintf  (’print  -dpng  Zs'/od.temperature  ’  ,project_name,k)) ; 
k  =  k+1; 

end 
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Pressure  (kPa)  at  time  t  =  2  ms 


Pressure  (kPa)  at  time  t  =  4  ms 


(i)  t  =  2  ms. 


(ii)  t  =  4  ms. 


Pressure  (kPa)  at  time  t  =  6  ms 


Pressure  (kPa)  at  time  t  =  8  ms 


(iii)  t  =  6  ms. 


(iv)  t  =  8  ms. 


(v)  t  =  10  ms. 


(vi)  t  =  12  ms. 


Figure  4-  Two-dimensional  blast  propagation;  pressure  field  surfaces. 
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(i)  t  =  2  ms. 


(ii)  t  =  4  ms. 


Temperature  (K)  at  time  t  =  6  ms 


Temperature  (K)  at  time  t  =  8  ms 


0  0  0  0 


(iii)  t  =  6  ms. 


(iv)  t  =  8  ms. 


Temperature  (K)  at  timet  =  10  ms 


Temperature  (K)  at  timet  =  12  ms 


(v)  i  =  10  ms. 


(vi)  t  =  12  ms. 


Figure  5:  Two-dimensional  blast  propagation;  temperature  field  surfaces. 
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6.4  Three-dimensional  blast  propagation 

Consider  three-dimensional  explosion  in  air  occupying  a  block  with  dimensions  25 
by  10  by  20  m  with  the  initial  conditions  similar  to  the  previous  section.  The  charge  is 
positioned  at  a  corner  of  the  block.  Both  reflective  and  non-reflective  boundary  conditions 
are  imposed  on  the  walls  of  the  block  to  model  blast  propagation  in  a  street  canyon  under 
the  assumption  that  the  charge  is  positioned  at  the  ground  level  at  the  centre  line  of  the 
street. 

The  following  MATLAB  script  simulates  this  problem  on  a  uniform  mesh.  The  pressure 
and  temperature  fields  for  increasing  times  are  presented  in  Figures  6  and  7.  It  is  clearly 
seen  the  formation  of  a  secondary  shock  wave. 


"/X  3D  blast  propagation  in  a  street  canyon 
project_naine  =  ’blastSd’; 

dispCsprintf  ( ^Project  name:  7«s\project_name)) ; 
gds  =  gas_dynamics_solver ; 

R  =  287.0;  KelvinO  =  273.0; 

TO  =  KelvinO+20.0;  pO  =  l.OleS;  rO  =  pO/TO/R; 

T1  =  KelvinO+2000.0;  pi  =  ll*p0;  rl  =  pl/Tl/R; 

CFL  =  0.5;  save_time  =  0.005;  kmax  =  6; 

/s'/o  Pre-process  if  necessary 
k  =  0; 

if  "exist (sprintf ( ’ZsZd.mat ^ ,project_name ,k) , ’ f ile O 
'/o  Generate  mesh 

s.LX  =  25.0;  s.LY  =  10.0;  s.LZ  =  20.0; 
s.NX  =  100;  s.NY  =  40;  s.NZ  =  80; 

dx  =  s.LX/s.NX;  dy  =  s.LY/s.NY;  dz  =  s.LZ/s.NZ;  s.dx  =  min( [dx,dy ,dz] ) ; 
s.x  =  -dx/2:dx:s.LX+dx/2;  s.y  =  -dy/2 : dy : s .LY+dy/2;  s.z  =  -dz/2 : dz : s . LZ+dz/2 ; 
NC  =  (s.NX+2)*(s.NY+2)*(s.NZ+2) ; 

NF  =  (s.NX+l)*s.NY*s.NZ+s.NX*(s.NY+l)*s.NZ+s.NX*s.NY*(s.NZ+l) ; 
dispCsprintf  (’“/od  facets,  “/od  cells  ’  ,NF,NC-8)  )  ; 

s.dv  =  (dx*dy*dz) *ones (NC , 1) ;  s.fc  =  zeros(NF,2);  s.ds  =  zeros(NF,3); 
s.BE  =  zeros  (s  .NY*s  .NZ ,  1)  ;  s.GE  =  zeros  (s  .NY*s  .NZ,  1)  ;  7,  east  BC 

s.BW  =  zeros  (s  .NY*s  .NZ ,  1)  ;  s.GW  =  zeros  (s  .NY*s  .NZ ,  1)  ;  7o  west  BC 

s.BS  =  zeros  (s  .NX*s  .NZ,  1)  ;  s.GS  =  zeros  (s  .NX*s  .NZ,  1)  ;  7o  south  BC 

s.BN  =  zeros  (s  .NX*s  .NZ ,  1)  ;  s.GN  =  zeros  (s  .NX*s  .NZ ,  1)  ;  7o  north  BC 

s.BB  =  zeros  (s  .NX*s  .NY,  1)  ;  s.GB  =  zeros  (s  .NX*s  .NY,  1)  ;  7o  bottom  BC 
s.BT  =  zeros  (s  .NX*s  .NY,  1)  ;  s.GT  =  zeros  (s  .NX*s  .NY,  1)  ;  7o  top  BC 
s.map  =  @(ix,iy,iz)  sub2ind( [s .NX+2 , s . NY+2, s .NZ+2] , ix , iy , iz) ; 
i  =  0;  7«  facet  index 
for  ix  =  l:s.NX+l 

for  iy  =  2:s.NY+l 

for  iz  =  2:s.NZ+l 

i  =  i+1;  s.ds(i,:)  =  [dy*dz,0,0]; 
s.fc(i,l)  =  s.map(ix,iy,iz) ; 
s.fc(i,2)  =  s.map(ix+l,iy,iz) ; 

end 

end 

end 

for  ix  =  2:s.NX+l 

for  iy  =  l:s.NY+l 

for  iz  =  2:s.NZ+l 

i  =  i+1;  s.ds(i,:)  =  [0,dx*dz,0]; 
s.fc(i,l)  =  s.map(ix,iy,iz) ; 
s.fc(i,2)  =  s.map(ix,iy+l,iz) ; 

end 

end 

end 

for  ix  =  2:s.NX+l 

for  iy  =  2:s.NY+l 
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for  iz  =  1 : s . NZ+1 

i  =  i+1;  s.ds(i,:)  =  [0,0,dx*dy] ; 
s.fc(i,l)  =  s.inap(ix,iy,iz) ; 
s.fc(i,2)  =  s.inap(ix,iy,iz+l) ; 

end 

end 

end 

for  iy  =  2:s.NY+l 

for  iz  =  2:s.NZ+l 

i  =  sub2ind( [s .NY, s .NZ] , iy-1 , iz“l) ; 

s.BE(i)  =  s .map(2, iy , iz) ;  s.GE(i)  =  s .mapCl , iy , iz) ; 

s.BW(i)  =  s.map(s.NX+l,iy,iz) ;  s.GW(i)  =  s .map(s . NX+2, iy , iz) ; 

end 

end 

for  ix  =  2:s.NX+l 

for  iz  =  2:s.NZ+l 

i  =  sub2ind( [s .NX , s .NZ] , ix-1 , iz-1) ; 

s.BS(i)  =  s.map(ix,2,iz) ;  s.GS(i)  =  s.inap(ix,l,iz) ; 

s.BN(i)  =  s.map(ix,s.NY+l,iz) ;  s.GN(i)  =  s .map(ix, s . NY+2, iz) ; 

end 

end 

for  ix  =  2:s.NX+l 

for  iy  =  2:s.NY+l 

i  =  sub2ind( [s .NX , s .NY] , ix-1 , iy-1) ; 

s.BB(i)  =  s.map(ix,iy,2) ;  s.GB(i)  =  s.inap(ix,iy,l) ; 

s.BT(i)  =  s.map(ix,iy,s.NZ+l) ;  s.GT(i)  =  s .map(ix, iy , s . NZ+2) ; 

end 

end 

%  Initial  conditions 

xa=0.0;  ya=0.0;  za=0.0;  a  =5.0;  b=  100.0; 
density  =  @(x,y,z) . . . 

0.5*((rl+r0)  +  (rl-r0)*tainh(.  .  . 
b*  (a-sqrt  (  (x-xa)  “2+(y-ya)  “2+(z-za)  ''2)  )  )  )  ; 
pressure  =  @(x,y,z)... 

0.5*((pl+p0)+(pl-p0)*tanh(. . . 
b*  (a-sqrt  (  (x-xa)  “2+(y-ya)  “2+(z-za)  ''2)  )  )  )  ; 
s.r  =  zeros(NC,l);  s.p  =  zeros(NC,l);  s.v  =  zeros(NC,3); 
for  ix  =  l:s.NX+2 

for  iy  =  l:s.NY+2 

for  iz  =  1 : s . NZ+2 

i  =  s .map(ix, iy , iz) ; 

s.r(i)  =  density(s.x(ix) ,s.y(iy) ,s.z(iz)) ; 
s.p(i)  =  pressure(s.x(ix) ,s.y(iy) ,s.z(iz)) ; 

end 

end 

end 

%  Specify  Riemann  solver 

s. gamma  =  1.4;  rs  =  riemann_solver(s .gamma) ; 

s. state  =  rs . interf ace_state ;  s.t  =  0.0;  s  =  gds . initialise (s) ; 
eval (sprintf  ( ’ save  7«s7,d  -struct  s’ ,project_name,k)) ; 

end 


7«7o  Find  the  latest  restart  file 

while  exist (sprintf(’7»s7od. mat’  ,project_name ,k) ,  ’file’) 
k  =  k+1; 

end 

k  =  k-1; 

7«7o  Load  start/restart  file 
s  =  load(sprintf  (’7os7«d’ ,project_name,k)) ; 
while  k  <  kmax 

elapsed_time  =  cputime; 
k  =  k+1; 

tf  =  s.t  +  save_time; 
while  s.t  <  tf 

s  =  gds .update_rates(s) ; 

dt  =  CFL*s.dx/max(s.ia+abs(s.ivn)) ; 
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dt  =  miii(dt,tf-s.t)  ;  s.t  =  s.t+dt; 
dispCsprintf  ( =  '/og’jS.t)); 
s.r  =  s.r+s.r_rate*dt; 
s.e  =  s . e+s . e_rate*dt ; 
s.u  =  s.u+s.u_rate*dt; 

7o  reflective  east  BC  (symmetry  plane) 
s.r(s.GE)  =  s.r(s.BE);  s.e(s.GE)  =  s.e(s.BE); 
s.u(s.GE,l)  =  “S.u(s.BE,l) ;  s.u(s.GE,2)  =  s.u(s.BE,2); 
s.u(s.GE,3)  =  s.u(s.BE,3); 

7  non-ref lective  west  BC 

s.r(s.GW)  =  s.r(s.BW);  s.e(s.GW)  =  s.e(s.BW); 
s.u(s.GW,l)  =  s.u(s.BW,l);  s.u(s.GW,2)  =  s.u(s.BW,2); 
s.u(s.GW,3)  =  s.u(s.BW,3); 

7  reflective  south  BC  (symmetry  plane) 
s.r(s.GS)  =  s.r(s.BS);  s.e(s.GS)  =  s.e(s.BS); 
s.u(s.GS,l)  =  s.u(s.BS,l);  s.u(s.GS,2)  =  -s .u(s . BS , 2) ; 
s.u(s.GS,3)  =  s.u(s.BS,3); 

7  reflective  north  BC 

s.r(s.GN)  =  s.r(s.BN);  s.e(s.GN)  =  s.e(s.BN); 
s.u(s.GN,l)  =  s.u(s.BN,l);  s.u(s.GN,2)  =  -s .u(s . BN, 2) ; 
s.u(s.GN,3)  =  s.u(s.BN,3); 

7  reflective  bottom  BC 

s.r(s.GB)  =  s.r(s.BB);  s.e(s.GB)  =  s.e(s.BB); 
s.u(s.GB,l)  =  s.u(s.BB,l);  s.u(s.GB,2)  =  s.u(s.BB,2); 
s.u(s.GB,3)  =  -s.u(s.BB,3) ; 

7  non-ref lective  top  BC 

s.r(s.GT)  =  s.r(s.BT);  s.e(s.GT)  =  s.e(s.BT); 
s.u(s.GT,l)  =  s.u(s.BT,l);  s.u(s.GT,2)  =  s.u(s.BT,2); 
s.u(s.GT,3)  =  s.u(s.BT,3); 

7  synchronise  physical  variables 
s  =  gds .update_states (s) ; 

end 

elapsed_time  =  cputime-elapsed_time ; 

disp(sprintf (’CPU  time:  7g  secs  /  7g  mins  /  7g  hours’,... 

elapsed_time ,elapsed_time/60 ,elapsed_time/3600) ) ; 
eval (sprintf ( ’ save  7s7d  -struct  s’ ,project_name,k)) ; 

end 

77  Post  process 

disp( ’Generating  PNG  images...’); 
set(clf , ’Renderer’ , ’ZBuffer’) ; 
k  =  0; 

while  exist(sprintf (’7s7d.mat’ ,project_name,k) , ’file’) 
s  =  load(sprintf (’7s7d’ ,project_name,k)) ; 
sx  =  s.x(s.NX+l);  sy  =  [s .y (2) , s .y (s .NY+1)] ;  sz  =  s.z(2); 
[x,y,z]  =  ndgrid(s.x,s.y,s.z) ;  w  =  zeros (size (x) ) ; 
for  ix  =  2:s.NX+l 

for  iy  =  2:s.NY+l 

for  iz  =  2:s.NZ+l 

i  =  s .map(ix, iy , iz) ; 
w(ix,iy,iz)  =  s.p(i)/1000; 

end 

end 

end 

slice( . . . 

y(2:end-l,2:end-l,2:end-l) , . . . 
x(2:end-l,2:end-l,2:end-l) , . . . 
z(2:end-l,2:end-l,2:end-l) , . . . 
w(2:end-l,2:end-l,2:end-l) , . . . 
sy,sx,sz) ; 

daspect ( [1 , 1 , 1] ) ;  axis  equal;  colorbar; 

title(sprintf ( ’Pressure  (kPa)  at  time  t  =  7g  ms ’ , 1000*s .t) ) 
eval (sprintf (’print  -dpng  7s7d_pressure’ ,project_name,k)) ; 
for  ix  =  2:s.NX+l 

for  iy  =  2:s.NY+l 

for  iz  =  2:s.NZ+l 

i  =  s .map(ix, iy , iz) ; 


28 


DSTO-TR~2139 


w(ix,iy,iz)  =  s .p(i) /s . r (i)/R; 

end 

end 

end 

slice( .  .  . 

y(2:end“l,2:end-l,2:end-l)  ,  .  .  . 
x(2:end-l,2:end“l,2:end-l) , . . . 
z(2:end-l,2:end“l,2:end-l) , . . . 
w (2 : end- 1, 2: end- 1,2: end- 1) , . . . 
sy,sx,sz) ; 

daspect ( [1 , 1 , 1] ) ;  axis  equal;  colorbar; 

title(sprintf  ( ’Temperature  (K)  at  time  t  =  “/g  ms’ ,1000*s.t)) ; 
evaKsprintf  (’print  -dpng  VoS'/od.temperature ’  ,project_name,k)) ; 
k  =  k+1; 

end 
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(i)  t  =  5  ms. 


(ii)  t  =  10  ms. 


(iii)  i  =  15  ms. 


(iv)  t  =  20  ms. 


(v)  t  =  25  ms. 


(vi)  i  =  30  ms. 


Figure  6:  Three-dimensional  blast  propagation;  pressure  field  slices. 
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(i)  t  =  5  ms. 


(ii)  t  =  10  ms. 


(iii)  i  =  15  ms. 


(iv)  t  =  20  ms. 


(v)  t  =  25  ms. 


(vi)  i  =  30  ms. 


Figure  7;  Three-dimensional  blast  propagation;  temperature  field  slices. 
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7  Conclusion 

The  two  reusable  M-files  constituting  a  Gas  Dynamics  Toolbox  for  MATLAB  are  suffi¬ 
cient  for  simulation  of  multi-dimensional  problems  for  blast  propagation.  The  comparison 
with  the  benchmark  solution  and  the  results  of  numerical  simulation  demonstrate  that  the 
code  provides  a  reasonable  solution  which  can  be  used  for  estimation  of  the  peak  pressure 
and  impulse  of  shock  waves. 

It  is  worthwhile  noting  that  the  described  numerical  scheme  is  an  implementation  of 
the  Finite  Volume  Method  which,  in  turn,  is  a  particular  realisation  of  the  Finite  Ele¬ 
ment  Method  when  both  the  shape  and  weighting  functions  coincide  with  the  character¬ 
istic  functions  of  cells  Wj.  The  Finite  Volume  Method  provides  a  conservative  numerical 
scheme  whose  solution  exhibits  highly  desirable  properties  of  conservation  of  total  mass, 
momentum  and  energy.  Though  the  described  first-order  accurate  Godunov-type  solver 
is  susceptible  to  numerical  diffusion  on  coarse  meshes,  the  total  quantities  such  as  the 
total  load  on  part  of  a  wall  (e.g.  windows  and  doors  of  a  building)  can  be  estimated  quite 
reliably  due  to  the  conservativeness  of  the  scheme. 
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